Paired phases and Bose-Einstein condensation of spin-one bosons with attractive 

interaction 
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We analyze paired phases of cold bosonic atoms with the hyper spin 5=1 and with an attractive 
interaction. We derive mean-field self-consistent equations for the matrix order parameter describing 
such paired bosons on an optical lattice. The possible solutions are classified according to their 
symmetries. In particular, we find that the self-consistent equations for the SO (3) symmetric phase 
are of the same form as those for the scalar bosons with the attractive interaction. This singlet 
phase may exhibit either the BCS type pairing instability (BCS phase) or the BEC quasiparticle 
condensation together with the BCS type pairing (BEC phase) for an arbitrary attraction Uo in the 
singlet channel of the two body interaction. We show that both condensate phases become stable 
if a repulsion U2 in the quintet channel is above a critical value, which depends on Uo and other 
thermodynamic parameters. 
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I. INTRODUCTION 

Superconductivity in metals is a consequence of pairing 
between electrons |l[ and formation of a new macroscopic 
coherent state made of electron pairs @. The micro- 
scopic Bardecn-Cooper-Schrieffer (BCS) theory of super- 
conductivity explains all key experimental features of the 
superconducting state. Recent developments in trapping, 
cooling, controlling, and detecting of atoms allowed to 
investigate superfluidity in neutral fermionic systems, cf. 
reviews: [!, Hjj. In particular, a crossover between Bose- 
Einstein condensation (BCS) limit, where fermionic pairs 
overlap significantly, and the BEC limit, where tightly 
bound pairs form a coherent state of Bose-Einstein con- 
densed bosons @, @ , was demonstrated experimentally 
0- Since both fermionic as well as bosonic atoms are 
available in experiments it is natural now to investigate 
pairing between bosons and a coherent superfluid state 
of paired bosons. 

Pairing and phase transitions of bosons with zero 
spin and with an attractive interaction was discussed in 
Ref. 8] in the context of superfluid helium four. It was 
found that collective excitations of a coherent condensed 
state of paired bosons can undergo another BEC type 
condensation into a one-particle condensed state, now 
known as the Evans-Rashid transition However, 
already in Ref. [Io[ it was found that both the homoge- 
neous coherent paired phase and the homogeneous phase 
due to the Evans-Rashid transition are unstable against 
a mechanical collapse. I.e., the bosons with an attractive 
interaction tend to form clusters of_particles. Moreover, 
extending the mean-field result of I 8HIOII by the leading- 
order fluctuation contributions [n], [l2[ or higher order 
corrections (l3j to the thermodynamic potential due to 
the interaction between particles does not stabilize ho- 
mogeneous phases of bosons with pairing potential. It 
is observed in [l4| that by confining bosons with paring 



interaction in a trap, which produces a gap between the 
ground and the excited states, one can protect the sys- 
tem from the mechanical instability. On the other hand, 
in Ref. [l5j a narrow region at finite temperatures on 
a phase diagram is found, where many-body effects can 
stabilize the homogeneous normal phase of bosons with 
pairing interaction. 

In the present paper we employ a mean-field theory 
to solve a problem of pairing between bosons with spin 
5=1 moving on an optical lattice. We show that the 
coherent BCS type phase of paired bosons induced by 
attractive interaction in the singlet channel is stable po- 
vided that the interaction in a quintet channel is repulsive 
and sufficienlty strong. Since bosons with the nonzero hy- 
perspin are available in experiments with cold atoms 0,0] 
and both sign and strength of the interaction between 
them can be tuned by the optical Feshbach resonance 
[l6| it is interesting to look for such a paired bosonic 
state in a laboratory. 

The paper is organized as follows: In Section II we in- 
troduce the model. In Section III spinor pair condensed 
phases are classified by the matrix BCS type order pa- 
rameter and the mean-field theory solution of the Hamil- 
tonian is presented. Section IV is devoted to the symme- 
try classification for the phases characterized by the com- 
plex matrix order parameter following the standard ap- 
proaches [l7j - [l9j . It should be noted that the finite spin of 
bosons with repulsive interaction in all channels [lo, Hl[ 
leads to many nontrivial grou nd and excited states of 
the spinor condensates [J, |221 . They include topologi- 
cally nontrivial phases pjl, H3 1 skyrmion excitations 
[25(, and even nonabelian vortices [26|. The symmetry 
classification presented in this Section may help in the 
future for detailed analysis of non-trivial excitations in 
paired phases of bosons. In Section V we present numer- 
ical solutions to the mean-field equations and in Section 
VI stability of particular phases is discussed. Conclu- 
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sions are in Section VII and details on our derivations 
are presented in Appendices. 

II. THE MODEL 

We consider the Hubbard model for spin-one bosons, 
which are trapped on an optical lattice. The grand 
canonical Hamiltonian [2?], 

H = H a + H int - fiN (2.1) 

contains the kinetic part H° and the interaction part 
H mt . We introduce the chemical potential /i to fix the 
average total number of bosons in the lattice. The kinetic 
part is 

H " = ~t E 6 U" ( 2 - 2 ) 

where bi a (^L) is an annihilation (creation) operator of 
a boson at the lattice site i with spin a — —1,0, or 1, 
and denotes the summation over nearest neighbor 
sites. We also introduce here the hopping integral t. We 
absorb a constant single site occupation energy into the 
definition of /x. The effective parameter of our model t 
can be derived from the microscopic details of the opti- 
cal lattice 0, , assuming that the lattice site orbitals 
correspond to localized Wannier functions with one level 
per site. We neglect the harmonic trap confinement in 
the following. 

The interaction part H lnt is constructed under an as- 
sumption that the total spin of the system is conserved 
and that the interaction amplitude is local [13, [13] ■ While 
the first requirement is natural due to general conser- 
vation laws, the second assumption is justified for cold 
atoms due to their neutrality and short-range character 
of interacting forces. The total spin S of two interacting 
spin-one bosons attains three possible values S — 0, 1, 2. 
Because of the bosonic symmetry of the wave functions, 
in the presence of the local interaction only 5 = and 
5 = 2 terms contribute in H mt . The resulting inter- 
action amplitudes Us are proportional to the scattering 
lengths as for each S channel [3]. Following Ref. [3l[ we 
write H int in terms of the number operator m and spin 
operator Si at the lattice site i: 

H mt = |E " 1) + | E( S ' " 2 «' (2-3) 

i i 

were g n = (2U 2 + U )/3 and g s = (U 2 - U )/3. Both 
the hopping integral t and the interaction strength Us 
can be tuned to become of comparable magnitude by 
manipulating the laser light producing the optical lattice. 

In this paper we are interested in the effects of at- 
tractive interaction giving rise to pairing between spin- 
one bosons. We introduce an auxiliary annihilation (and 
creation) operator for a Cooper pair of bosons B^ — 
J2oo' C^^fbiabjcr' , where C^'Jf is the Clebsch-Gordan 



coefficient for the total spin S and with spin projection 
M . The explicit form of the pair operators can be found 
in Ref. [22j ■ The Hamiltonian (|2.3[) takes a new, compact 
form 

BP* = J2 (u B^B°f + U 2 j2 Bl^BlA , 

i \ M=-2 / 

(2.4) 

which is more appropriate here since it shows directly 
all structures of bosonic pair correlations and hints to 
possible order parameters. In the next Section we solve 
the model (|2~Tj) with (j2~2j) and (|2~l3"|) within a Hartree- 
Fock mean-field approximation (MFA) and discuss pos- 
sible condensed phases of bosonic Cooper pairs. 

III. MEAN-FIELD APPROXIMATION 

The kinetic part (|2.2p of our model Hamiltonian is di- 
agonal in the momentum representation 

H° =J2^Lb k<7 , (3.1) 

ka 

where b£ (bka) is the creation (annihilation) operator 
for a particle with the lattice momentum k and the sin- 
gle particle kinetic energy is denoted by We keep a 
general form of the dispersion relation in our deriva- 
tion of the self-consistent equations and use a specific 
model later in Section |Vj 

The construction of the appropriate mean-field Hamil- 
tonian can be done within a textbook rule [32j by split- 
ting two-body operators into paring operators and their 
non-vanishing expectation values. This approximation 
consequently neglects fluctuations. The form of our 
model Hamiltonian (|2.4[) suggests the following choice for 
the pair expectation value 

A S,M = {B S,M h (32) 

which can also be expressed as A S ' M = Ylncr' C„'^f A CT<7 ' , 
with A aa > = -^-J^ki^kcrb-ka')- The expectation val- 
ues are taken at thermal equilibrium with the inverse 
temperature (3 and N s denotes the number of lattice 
sites. We also allow for nonzero normal density expecta- 
tion values by defining an average site occupation matrix 
n aa > = -i^Yskibk&bka')- Throughout the paper we deal 
with quantities described by 3 x 3 matrices in the spin 
index, such as n aa i or A ff(T < . Therefore, we introduce 
here a more compact matrix notation h and A for those 
quantities. 

Within MFA |], [Hi the Hamiltonian (JH3J) takes the 
following form 

#mf = E ( 6 L«W^' + WkA°°'bl k *> +h.c.))-E N s . 

k<j<j' 

(3-3) 
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In the above equation a matrix valued order parameter 
appears 

2 

A = U C 0fl A a '° + U 2 &' M A 2 > M , (3.4) 

M=-2 

which describes spontaneous symmetry breaking due to 
BCS-type paring. The quantity 



= 2[u Q C ' nC ' a + U 2 C 2M nC 2 > M 



(3.5) 



M=-2 



describes the effective Hartee-Fock potential. Note that 
the Clebsch coefficients for fixed 5, M are also repre- 
sented by a 3 x 3 matrix C S ' M . We keep the additive 
constant 
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2Tr(n T C ^ n(7^ ) + |Tr(A(7 ' )| 2 
Tr(n 2 + A f A) + n 2 



(3.6) 



which is necessary in the discussion of phase stabilities 
presented in the Section IVll 

We finally arrive at the mean field self-consistent equa- 
tions by calculating the normal n aa > and anomalous 
A CTCr ' averages in the grand canonical ensemble with the 
quadratic interaction Hamiltonian (|3.3I) . This procedure 
is equivalent [32| to the requirement of attaining a mini- 
mum of the free energy with the Hamiltonian (13.31) , when 
A and w are variational parameters. The technical de- 
tails of the derivation are given in the Appendix [X] Here 
we present the final result obtained from (|A3j) and (|A9[) 



t + h* 
A* 



-A 
— n 



where Mk is a Bogoliubov-de Gennes matrix 

_ ( :- M)l + w A 

-A* -{^-^t-w* 



(3.7) 



(3.8) 



For the purpose of solving the self-consistency equa- 
tions (|3.7p and (|3.8p in practice it is convenient to sim- 
plify the expression for A given in (|3.4p and for w in 
(|3.5[) . With the help of general algebraic identities 33 1 
applied to the matrices A and n we arrive at 

A = (U - U 2 )C ' Tt(C '°A) + U 2 A, (3.9) 

w = 2(U - U 2 )C°>°nC ' + U 2 {n T + nl), (3.10) 
where all C S,M have been eliminated except of C°'°. 



IV. SYMMETRY CLASSIFICATION OF 
ORDERED STATES 

The accepted strategy, which allows to classify the 
solutions for the matrix order parameter from the self- 
consistent equations, relays on symmetry considerations 



jl8| . The symmetry arguments alone allow to identify 
stationary states of the free energy, as it was done re- 
cently for the spinor condensates [JjJ HH . Here we need 
not only to identify the symmetry classified states, but 
we also want to investigate the phase diagram as a func- 
tion of the interaction parameters. Therefore, we have to 
compare free energy of symmetry classified phases to find 
the minimal one. In the investigation of superfluid 3 He it 
was observed [l?} , but not strictly proven, that the phase 
possessing the highest remaining symmetry corresponds 
indeed to a local, and very often to the global free energy 
minimum. 

We follow the standard symmetry classification ap- 
proach. We start by determining the highest allowed 
symmetry phase, and then we consider the solutions with 
a lower symmetry. For the sake of completeness of the 
presentation we give below a more detailed account of 
this derivation. We will use the classification introduced 
in this Section to determine numerically the phase di- 
agram, by solving the non-linear mean field equations 
within a given symmetry class. 

The full symmetry of our system (in a generic case 
U 2 7^ Uq) involves the gauge and the spin rotation sym- 
metry, so it is U(l) x 5*0(3). This symmetry is smaller 
then in the superfluid 3 He case, which has U(l) x 5*0(3) x 
50(3) symmetry group. The possibility of breaking the 
gauge invariance is crucial in our search of the phases 
with pairing. Symmetry of our system allows the gauge 
symmetry to be broken not only independently, but also 
in a combination with the spin symmetry operation. 
Thus we have to consider also the possibility of gauge- 
spin symmetry breaking, similar to superfluid 3 He. 



A. Symmetry transformations 

We start the discussion of symmetry with the global 
U(l) gauge symmetry transformation bk„ — > e % ^bk a , 
where %j) is a constant phase. The single site occupation 
matrix h is gauge invariant, so from (|3.5p it follows that 
w is gauge invariant as well. The pair expectation value 



transforms as A 



,2i 



^A, which substituted to ([3^3]) 
leads to the order parameter transformation A — > e 2l ^A. 
It is easy to check that this gauge transformation is a 
symmetry of our mean field equation p. 71) with (I3.8p . 

The spin rotation SO(3) is described by a unitary ma- 
trix f, which acts as follows: bka ~> So-' r a-a'bka'- The 
general rotation matrix f can be parameterized by three 
Euler angles of elementary rotations generated by three 
components of the spin-one operator. From the defini- 
tions of the averages h and A we obtain the transforma- 
tion rules 



n — > r nr 



A 



? Af T . 



(4.1) 



The above transformations substituted to (|3 . 10[) and 
(I3.9P give the following spin rotation of the effective po- 
tential and the pairing order parameter: 

A^^fAr)^, (4.2) 



w 



fwf\ 
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where r) — ^/3C°'°. We have used the identity 
f]f*f] — f, which follows from the explicit form •q aa > = 
— {—iy5 a - a i. One can check that the right-hand side 
of the mean field equation Eq. (|3.7p consequently trans- 
forms as Mfc — > RMkR\ with a unitary R — diag(f , f *), 
where diag stands for a block diagonal matrix. The 
left-hand side of this equation transforms upon (|4.1[) in 
the same manner, thus verifying the SO(3) spin rotation 
symmetry of our mean-field formulation. 

B. Continuous symmetry phases 

No broken symmetry. The requirement of invariance 
upon the full symmetry transformation U(l) x SO(3) ap- 
plied to w and A gives as the only solution w — wl and 
A = 0, where w — ^(U + 5U 2 )n. The single site density 
is n and the occupation matrix reads h = ^nt. This 
describes a free boson gas with a renormalized chemi- 
cal potential due to the Hartree-Fock treatment of the 
contact interaction. 

Singlet phase. The highest possible symmetry phase 
with non-zero pairing amplitude arises when we break 
the U(l) gauge symmetry, but leave the spin rotation 
symmetry. We derive from the invariance condition 

Afj = rhfjf* (4.3) 

that for a general f the order parameter has to be 
Ar) = Al with some complex A and w = tol, with w 
the same as in the free case discussed above. Going back 
to Eq. (|3.4[) we find that the expectation values of the 
bosonic pair operators A S,M are non-zero only for 5 = 
in this SO (3) symmetric phase. We will call this phase 
the singlet phase as pairing happens only in the singlet 
channel, with the finite order parameter A = UqC '°A ' . 

Quintet phase. We search now for paired phases, 
which allow for non-vanishing S = 2 (i.e. quintet) com- 
ponents of the order parameter. The simplest way to 
achieve this is by lowering the spin rotation SO (3) sym- 
metry to an axial L^(l) symmetry. We choose an arbi- 
trary quantization axis and express the spin rotations 
around this axis as f (cp) = e i(pS * , with some angle <p and 
S z = C 2 - 2 — C 2 '~ 2 . We require now a more general spin- 
gauge invariance condition for the order parameter 

Af/ = e 2ljf, r(ip)Afjr(tp)^ (4.4) 

where the gauge symmetry breaking phase tp can now 
depend on the spin rotation angle ip. We obtain three 
different solutions, which are presented below: 

U(l) s ,- V : i/) = -<p A = AC 2 ' 2 , (4.5a) 
: *l> = -<p/2 A = AC 2 ' 1 , (4.5b) 
U{l)s, : V = A = AC* 0,0 + A'C 2 ' . (4.5c) 

We follow the notation of Rcf. 17] to label the above 
spin-rotation breaking axial phases. Remaining solutions 



with +<p, and +p/2 can be obtained by changing the 
direction of the quantization axis, so they do not describe 
a different symmetry phase. 

The only non-zero pair expectation amplitude is A 2,2 
for the U{l)s z - V phase and A 2,1 for U(l)s,-^, which 
follows from the comparison of A definition in (|3.4[) with 
the result (|4.5j) . In these two axial phases the symme- 
try allows for pairing only in the quintet channel. The 
remaining U{l)s z phase has a mixed singlet-quintet pair- 
ing order parameter, which has to be parameterized by 
two (complex) numbers A and A'. 

The Hartree-Fock potential w in all the axial phases is 
restricted by the symmetry to be diagonal. This brings a 
possibility of magnetic order, coexisting with the pairing, 
marked by spin rotation symmetry breaking in the spin 
dependent site occupation. 



C. Discrete symmetry phase 

Within only 3x3 matrix representations one cannot 
construct the icosahedral or octahedral symmetry, with- 
out allowing for generation of all possible rotations. The 
biggest non-trivial discrete symmetry is thus T - the 
symmetry group of tetrahedron without reflections. The 
set of group generators can be explicitly expressed as 
{t,e l ^~ s ' )e i^(S„+V2S x )/V3)^ wnere we use 3 x 3 matrix 
representation of spin one with S x = C 2,1 + C 2 '^ 1 . Sub- 
stituting these generators for f in the invariance condi- 
tion (|4.4I) we obtain as the only solution ip = and 

A = A(C 2 - 2 + V2C 2 -- 1 ). 



V. MEAN-FIELD SOLUTION FOR SINGLET 
PHASE 

The singlet phase, introduced from the symmetry argu- 
ments in Section flVBl is our natural candidate for a phys- 
ically attainable phase. The system in the singlet phase 
has a maximal remaining symmetry of all the phases with 
non-zero pairing. The singlet phase is unitary, meaning 
that the matrix order parameter A is proportional to a 
unitary matrix. Stable phases of liquid 3 He were previ- 
ously found to be unitary [l7j as well. 

Simple form of the order parameter A = Afj in the sin- 
glet phase leads to an identity M% = e 2 l6x6, where the 
Bogoliubov-de Gennes matrix Mj, was defined in (13. SJ) , 
The quasi-particle excitation energy 

e k = - A* + wf - |A|2 (5.1) 

is a triple degenerate eigenvalue of as defined in Eq. 
(IA6|) . We can now directly calculate the generalized oc- 
cupation factor in the mean field equation (|3.7|) 

(l-e-W*)- 1 = f(e k )M k + ±, (5.2) 
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with /(efc) = 2 &k k ' ^ e recau that depends on 

w and A, which are related to h and A: 



w = j(U 
A = U A, 



+ 5U 2 )n, 



(5.3a) 
(5.3b) 



in the singlet phase, as obtained in section IIVBI We 
substitute (|5.3[) into in (|5.2I) and then equate to the 
left-hand side of Eq. (I3.7p . The resulting self-consistent 
equations in the singlet phase take a simple form 



n 
3 

1 



(5.4a) 
(5.4b) 



The first equation provides a relation between the av- 
erage occupation n and the chemical potential fj,, while 
the second guarantees a non-zero paring amplitude. The 
form of this second equation is similar to the gap equation 
in the BCS theory, but with a different function /(e/c) due 
to boson statistics of condensating quasiparticles. 

Interestingly, these equations are formally equivalent 
to the one obtained in the case of scalar attracting bosons 
in Ref. [l(| • The only difference is that the optical lattice 
provides a natural ultraviolet cutoff in our model. 

The existence of the BCS type singlet solutions de- 
pends only on the strength Uq of attraction in the singlet 
channel and is insensitive to scattering in the quintet 
channel. We will show in the next Section that the sin- 
glet paired phase of attracting bosons can be stabilized 
by a repulsive quintet interaction. This is in a marked 
difference to the scalar case, where the system always 
undergoes a mechanical collapse before reaching Evans- 
Rashid transition (Toj . 

We note that the quasiparticles in BCS type bosonic 
condensate may undergo a statistical (Bose-Einstein) 
condensation [H, [lj|. The transition occurs when the 
excitation spectrum in Eq. (|5.ip becomes gapless [35| . 
The singular condition ek=o = can be satisfied in a 
thermodynamic limit for 



|A| = &= 



fj. + w, 



(5.5) 



which fixes the chemical potential similarly to a standard 
BEC. We separate the k = terms to obtain 



/(<**)= 



"BEC 

3IAI 



(5.6) 



where we introduce tt-bec ~~ a finite average density for 
quasi-particles with k — only. The self-consistent equa- 
tions for the BEC quasiparticle phase follow 



n _ TlBEC 
3 ~~ 3 



fc/0 



+ |A| 2 /(e fc ) 



1 _ "-BEC 

U ~ 



3IAI 



fe^O 



(5.7a) 
(5.7b) 




FIG. 1: The diagram of the singlet phase presented in the den- 
sity n - temperature T coordinates at Uo = —0.33W. The 
solid line denotes the boarder of BCS type phase with pair- 
ing ('NP' stands for no-pairing, normal phase), the dashed 
line marks the borderline of BEC quasiparticle condensate. 
For {72/1 f/o | = 0.59 both phases in (c) (blank) region ful- 
fill the standard thermodynamic stability conditions, but 
are unstable in (a) (grey) and (b) (light-grey) regions. For 
U2/\Uo\ = 0.64 the regions (b) and (c) are thermodynami- 
cally stable, (a) region is unstable. 




FIG. 2: The diagram of the singlet phase for strong interac- 
tion Uo = —W in the density n - temperature T coordinates. 
The zero temperature critical density is n* w 0.15. The solid 
line denotes the border of BCS type paired phase, the dashed 
line marks the borderline of BEC quasiparticle condensate. 
Stability regions (a), (b) and (c) are defined in the same way 
as in Fig. [T] 



when we substitute the decomposition (|5.6[) into (157 
The chemical potential fi is fixed by (|5.5p . so hbec be- 
comes a new thermodynamic parameter, which we have 
to determine. Therefore, we distinguish two different 
phases: i) BCS phase where nsEC = and A 7^ 0, and ii) 
BEC phase where tibec 7^ and A / 0. The condition 
for the BCS/BEC borderline is obviously u-qec — > 0, it is 
when (|5.7p reduces to (|5.4[) . We note that the transition 
between BCS and BEC in the boson case cannot be inter- 
preted as being a counterpart of the BCS-BEC crossover 
known in the fermionic condensed systems 0, Q . 
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In order to solve numerically Eqs. (|5.4I) and (15. 7|) 
we only need to provide the density of states p(E) = 
^2i~8(E — in the optical lattice. We assume 
a simple elliptic model for the density p{E) = 
^VW 2 ) 2 - E 2 , where W = (32tt) 2 ^H is chosen to 
fit a low-energy density profile obtained from the disper- 
sion relation 

We are able to gain some analytical insight into the 
solution for this BEC singlet phase. The integrals ap- 
pearing in (|5.7[) can be performed in the limit T — > 
leading to 



n = riBEC + ~ (1 
W _ 2rt B EC 4 



ui 2 ) arctan -h= 



(1 + lu) arctan A= 



(5.8a) 
, (5.8b) 



where u = y^-. These two non-linear algebraic equa- 
tions determine tibec and |A|/W for a given n. We de- 
fine a critical density n* by setting tibec = in (|5.8p . 
The BEC condensate solution with finite tibec exists for 
densities larger than n* at T = 0. For weak interactions 
|i/o|/W < \ we find only n* = solution, which means 
that there is only the BEC phase at T = 0. For stronger 
interactions |?7o|/W > \ we find a region of BCS phase 
extending down to zero temperature. We present the 
resulting finite temperature phase diagram in the weak 
interaction case in Fig. [T] for a fixed attractive interac- 
tion Uq = — Q.33W. The situation with strong interac- 
tion is illustrated in Fig. [2] for Uq — —W. Additionally, 
one can obtain an analytic solution n* — -^\Uo\/W for 
\U \/W>1. 



VI. STABILITY 

We discuss below the standard thermodynamic stabil- 
ity conditions expressed by: i) positivity of pressure p, 
ii) positivity of constant volume specific heat cy , and iii) 
positivity of isothermic compressibility kt (for the cal- 
culation see Appendix B). 

Singlet phase. We find that cy is positive in the singlet 
phase and is independent of U 2 . The pressure p and 
the inverse compressibility k^ 1 have a following linear 
dependence on U%: 



p(U ,U 2 )=p(Uo) + U 2 ^, 



10wf 
9a 3 



(6.1a) 
(6.1b) 



This means that for any point on the phase diagram in 
Fig. Q] or Fig. [2] we can find U 2 large enough to stabilize 
the BCS or BEC singlet phase. The shadowed regions in 
these figures exemplify stability for \U 2 \/Uo = 0.59 and 
0.64, respectively. 

We illustrate the nature of subsequent transition by 
plotting the specific heat Cy in Fig. [3] We show the 
specific heat dependence on the temperature T at fixed 



Uo/W=-1.0 

-0.33— - 
BEC/BCS □ 
BCS/NP a 



3 
2 : 




0.4 0.8 1.2 

k B T/W 



FIG. 3: The specific heat versus temperature for a fixed den- 
sity n — 1.5. The solid line illustrates the sequence of tran- 
sitions for interaction strength Uo = —W, the dashed line is 
for Uo = -0.33W. 



density n = 1.5, which is representative both for strong 
and weak attraction and contains all three phases: non- 
pairing, BCS-like and the BEC quasiparticle condensate. 
The plot does not depend on the strength U 2 , provided 
it is strong enough to stabilize the phases. The specific 
heat exhibits a jump at the onset of pairing (marked by 
a triangle in Fig. [3]), indicating that the correspond- 
ing transition is of second order. The transition to BEC 
quasiparticle condensate is contiunuos with a cusp in the 
specific heat dependence (marked by a square), the be- 
haviour being known in the usual Bogoliubov theory of 
BEC condensation @. 

We inspect further the details of stability lines shown 
in Fig. [2j We find generically two stable phases separated 
by an unstable one for the system at fixed temperature. 
The system will then have a tendency to spontaneously 
separate into the dense and dilute phases. We make this 
statement qualitative by considering the thermodynamic 
spinodal decomposition [36j into the dense BEC and di- 
lute BCS or normal phase. The results are presented in 
Fig. 21 The spinodal stability lines are redrawn from 
Fig. [2J region (b). The line marked by triangles is given 
by the compressibility condition k^ 1 = 0, while the one 
marked by the squares is given by the pressure p = 0. 
We find that the region denoted by light gray shadow- 
ing corresponds to a metastable state, which undergoes 
the spinodal dcomposition. The regions denoted BCS or 
BEC above the solid binodal line are stable against such 
a thermal fluctuation. This result indicates that the ther- 
mal fluctuations around the mean field solution do not 
change qualitatively our phase diagram, they are of im- 
portance at the vicinity of the stability borderlines. The 
role of quantum fluctuations is left for future research. 

We find that the most unstable point of our diagram 
is located at the BCS/BEC borderline (marked by the 
thick dashed line in Fig. @|. We can solve the following 
conditions p > and dp/dn > at zero temperature on 
this corssover line, corresponding to the density n* (com- 
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FIG. 4: The phase diagram of the singlet phase with the 
same parameters as in Fig. [21 The binodal line is indicated 
by a thin solid curve. The system in both shadowed regions 
(d) within BCS phase and (e) within BEC phase undergoes a 
spinodal decomposition. The BCS phase above (d) is stable, 
same as BEC to the left of (e) region. 



pare Fig. [2]). We thus find JTJ - the critical strength of 
repulsion in the quintet channel at which the whole BCS 
and BEC phases become stable. In the weak interaction 
regime 

which is valid for |Z7o|/W < \, while in the strong inter- 
action regime 

no = M(i+ 3arctan(l/^) \ 
2 20 V y/u - warctan(l/V^)y ' 

valid for |Z7q|/W > | and with oj calculated from 
Eq. (|5.8bj) at tibec = 0. For the interaction strength 
U2 > U% there is always a non-collapsing phase, which 
can be either normal, BCS or BEC homogeneous, or an 
inhomogeneous mixture of the dilute and dense phases. 

Other symetry phases. We have carried out a de- 
tailed, both analytical and numerical study of stability 
37] for the other phases. We find that both magnetic 
phases U(l)s z - V and U(l)s z ~ v /2 are mechanically un- 
stable. The tetrahedral T-phase with the quintet pairing 
occurs for U2 < and Uq > 0. It can be made thermody- 
namically stable by increasing the singlet repulsion Uq. 
We find however, that the ferromagnetic phase U(l)s z ~ v 
has lower free energy in the parameter regions, where 
T-phase becomes stable. Morover, we find that an in- 
finitezimal e > distortion of T-phase order paremctcr 
by a magnetic contribution A((l + e)C 2 ' 2 + V^C 2 ^ 1 ) 
leads to a lower free energy. We conclude that the T- 
phase is not even metastable, as it always corresponds to 
a saddle point of free energy. 



VII. CONCLUSIONS AND OUTLOOK 

In summary, we applied the Hartree-Fock mean-field 
approximation to solve a problem of pairing between 
bosons with spin S = 1 moving on optical lattices. The 
order parameter describing such paired bosons has a ma- 
trix form. Detailed classification of possible solutions ac- 
cording to their symmetries was presented. In particular, 
we found that the self-consistent equations for the SO(3) 
symmetric phase have the same form as those for the 
scalar bosons. We showed that the coherent BCS type 
phase of paired bosons induced by attractive interaction 
in the singlet channel is stable provided that the interac- 
tion in a quintet channel is repulsive. This finding might 
be usefull in experiments to stabilize bosons with attrac- 
tive interaction against mechanical collapse. 

The analyzed problem might be extended in the fu- 
ture in different ways. For example, it would be inter- 
esting to include local quantum correlation beyond the 
static Hartree-Fock approximation by using the bosonic 
dynamical mean-field theory developed recently [39j . An- 
other line of research is to investigate inhomogeneous ex- 
cited states of bosons with S = 1 in the BCS or BEC 
phases, i.e. there should be generalized vortex states in 
the condensed phases because of the high remaining sym- 
metry in the system. In the boson gas with hyper spin 
S = 2, where a very reach variety of spinor BEC for re- 
pulsive interactions have been proposed [401 ], we expect 
stabilization of at least some of many symmetry allowed 
phases for attractive interactions. 
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Appendix A: Diagonalization of the mean field 
Hamiltonian 

In this section we diagonalize our mean field Hamil- 
tonian (|3.3I) and derive Eq. (|3.7I) . We use a convenient 
notation [38j j , known as the Nambu notation in the stan- 
dard theory of superconductivity for fermions [32j 

*fc = (b k i, b k0 , b k -i, bl kv bl kQ , bl_ 1 ) T , (Al) 

where the superscript T denotes transposition. The 
bosonic canonical commutation relations rewritten with 
the Nambu spinor are (*t)/s] = i^z)aph,k> 

where a, (3 = 1, . . . , 6 and E z = diag(l, —1). Here E z is 
a diagonal matrix, the symbol l denotes the unity 3x3 
matrix. We express the model Hamiltonian in the mean 
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field approximation (|3.3[) as follows 

Hmf = H°+H™-iiN=±J2 *lS z A4* fc + const.. 

k 

The additive constant does not enter the calculation pre- 
sented in this Appendix. We introduce here the 6x6 
matrix M k resulting from the commutation [4>fe, Hmf] = 
M k & k - The explicit form of M k is given in Eq. (|3 . 8|) . 

With the Nambu spinor we write a compact expression 
for all normal and anomalous averages 

*EM> -ft.** J). <«> 

where the single site occupation n and the amplitude of 
Cooper pair condensate A were defined in Section Mil 
Our aim is now to compute the l.h.s. of the above equa- 
tion. This is easily done with a suitable Bogoliubov 
transformation performed on the mean field Hamiltonian 
(|3.3[) . We follow this route by introducing a new Nambu 
spinor T for quasiparticle excitations 

(A4) 

where u k and Vk contain coefficients to be determined 
below. We require the new spinor T to describe proper 
quasiparticles, so \T k ,HMF] = E k T k , where E k = 
diag(e fe ,-e fc ) and {e k )„„> = e ka 5 aa >. The eigenvalues 
efco- correspond to the excitation energy of quasiparti- 
cles in the BCS condensate. The components of T k 
have to fulfill the bosonic commutation relations, namely 
[(Tk) a , ( r l')/3] = (S z ) a ^5fc,fc'. The corresponding re- 
quirement for the coefficients of the Bogoliubov trans- 
formation (|A4[) leads to 

ft ?0 _1 =E,fji ^Ve,. (AS) 
V v fc «*/ \ v k u kJ 

Finally, we write the eigenvalue equation for the 6x6 
matrix M k 







( u k 


Vk\ 


= E k 













It follows from the hermicity of Y> z M k and the transfor- 
mation constrain (|A5|) that the eigenvalues of M k have to 
be real. Note also that the additional matrix E 2 enters 
our derivation due to the bosonic commutation relation, 
but is absent in the standard BCS formulation for the 
fermions. 

It turns out that we don't need to calculate explicitly 
u k and v k from Eq. (|A6|) as long as we are only inter- 
ested in the thermodynamic averages. The quasiparticle 
averages are particularly simple 

(T k T\,) S 2 - (1 - e"^*)- 1 ^. (A7) 



We transform this equation back to the original spinor 
<I> fc with the help of Eq. (|A4|I and we get 

(A8) 

which simplifies to 

(# fc *t)s, = (l-e-^) _1 - (A9) 

The above equation together with Eq. (|A3[> gives the final 
result of this Appendix. 



Appendix B: Calculation of thermodynamic 
parameters 

We calculate the grand canonical potential (per site) 

Q = -^lnTre~ mMF (Bl) 

by taking the trace Tr over all many-particle states of 
second-quantized Hmf as defined in (|A2|) . The general 
expression for the entropy per site S = — j£p\n is then 

S = ^(lnTre-^^ + (3 (H MF ) MF ) (B2) 

where (. . .) MF denotes the thermodynamic average with 
our model mean field hamiltonian. Using the results of 
diagonalization of bilinear Hmf derived in Appendix A 
we get 

S=^(-Hl-e-^) + J^j), (B3) 

while the constant introduced in (|3.3j) does not enter this 
expression. With the above explicit of S we write the 
final compact formula for fi 

n = -TS+-fc^(Z k -n)n k +Eo, (B4) 

k 

where n k — \J e k + |A| 2 /( e fc) — i is average occupation 
of a quasi-momentum state k. The constant has been 
now properly recovered and is entirely included in Eq as 
defined in (13.61) . For the singlet phase there is a simple 
scalar expression E = + 5U2 g U ° n 2 . With the grand 
canonical potential f2 expressed entirely in terms of the 
order parameter one calculates the pressure p — — Sl/a 3 , 
the specific heat per site cy = T^\ n and the inverse 
compressibility k^, 1 = n^-\r- 



9 



[1] L. N. Cooper, Phys. Rev 104, 1189 (1956). 

[2] J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Phys. 

Rev. 108, 1175 (1957). 
[3] D. Jaksch and P. Zoller, Annals of Physics 315, 52 

(2005). 

[4] I. Bloch, J. Dalibard, W. Zwerger, Rev. Mod. Phys. 80, 
885 (2008). 

[5] A. J. Leggett, in Modern Trends in the Theory of Con- 
densed Matter eds. A. Pekalski and R. Przystawa, pp. 13- 
27 (Springer- Verlag, Berlin, 1980). 

[6] P. Nozieres and S. Schmitt-Rink, J. Low Temp. Phys. 59, 
195 (1985). 

[7] C.A. Regal, C. Ticknor, J. L. Bohn, and D. S. Jin, Nature 

424, 47 (2003). 
[8] W. A. B. Evans and Y. Imry, Nuovo Cimento 63B, 155 

(1969). 

[9] W. A. B. Evans and R.I.M.A. Rashid, J. Low Temp. 

Phys. 11, 93 (1973). 
[10] H.T.C. Stoof, Phys. Rev. A 49, 3824 (1994). 
[11] E. Mueller and G. Baym Phys. Rev. A 62, 053605 (2000). 
[12] G. S. Jeon, L. Yin, S. W. Rhee, and D. J. Thouless, Phys. 

Rev. A 66, 011603, (2002). 
[13] M. Mannel, K. Morawetz, P. Lipavsky, New J. Phys. 12, 

033013 (2010). 

[14] P. Zin, B. Oles, M. Trippenbach, and K. Sacha, Phys. 
Rev. A 78, 023620 (2008). 

[15] A. Koetsier, P. Massignan, R. A. Duine, H. T. C. Stoof, 
Phys. Rev. A 79, 063609 (2009). 

[16] G. Thalhammer, M. Theis, K. Winkler, R. Grimm, and 
J. H. Denschlag, Phys. Rev. A 71, 033403 (2005). 

[17] C. Bruder, D. Vollhardt, Phys. Rev. B 34, 131 (1986). 

[18] D. Vollhardt and P. Woelfle, The superfluid phases of He- 
lium 3 (Taylor and Francis, London, New York, Philadel- 
phia, 1990). 

[19] S.-K. Yip, Phys. Rev. A 75, 023625 (2007). 
[20] T. Ohmi and K. Machida, J. Phys. Soc. Jpn. 67, 1822 
(1998). 

[21] T.-L. Ho, Phys. Rev. Lett. 81, 742 (1998). 
[22] M. Ueda, Y. Kawaguchi, "Spinor Bose-Einstein conden- 
sates", preprint \arXiv:100L2 072 
[23] M. Koashi and M. Ueda, Phys. Rev. Lett. 84, 1066 



(2000). 

[24] M. Ueda and M. Koashi, Phys. Rev. A. 65, 063602 

(2002) . 

[25] U. A. Khawaja, H. Stoof, Nature 411, 918 (2001); U. A. 

Khawaja, H. Stoof, Phys. Rev. A 64, 043612 (2001). 
[26] G. W. Semenoff, F. Zhou, Phys. Rev. Lett. 98, 100401 

(2007) . 

[27] M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. 

S. Fisher, Phys. Rev. B 40, 546 (1989). 
[28] D. van Oosten, P. van der Straten, and H. T. C. Stoof, 

Phys. Rev. A 63, 053601 (2001); Phys. Rev. A 67, 033606 

(2003) . 

[29] G. Mazzarella, Phys. Rev. A 73, 013625 (2006). 

[30] C. J. Pethick and H. Smith, Bose-Einstein Condensa- 
tion in Dilute Gases (Cambridge University Press, Cam- 
bridge, UK, 2002). 

[31] Imambekov Phys. Rev. A 68, 063602 (2003). 

[32] H. Bruus and K. Flensberg, Many-body quantum theory 
in condensed matter physics: an introduction (Oxford 
University Press, New York, USA, 2004). 

[33] The algebraic identities J2s M C s ' M Tx{C s ' M A) T = A 
and Y.s,mC S ' M MC S ' m ) t = (Tri)l, follow from the 
completness of total spin \S,M) basis. They are a ma- 
trix form of the ortogonality condition for the Clebsch- 
Gordan coefficients. 

[34] T.-L. Ho, and S.-K. Yip, Phys. Rev. Lett. 82, 247 (1999). 

[35] J. Cao, Y. Jiang, and Y. Wang, Europhys. Lett. 79, 
30005 (2007). 

[36] P.M. Chaikin, T.C. Lubensky Principles of condensed 
matter physics (Cambridge University Press, Cambridge, 
UK, 1995). 

[37] G. Pelka, PhD thesis, University of Warsaw (2010), in 
preparation. 

[38] H.T.C. Stoof, D.B.M. Dickerscheid, and K. Gubbels, Ul- 

tracold Quantum Fields, (Springer, Netherlands 2008). 
[39] K. Byczuk and D. Vollhardt, Phys. Rev. B 77, 235106 

(2008) . 

[40] R.B. Diener and T.-L. Ho, Phys. Rev. Lett. 96, 190405 
(2006). 



